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We describe the first axisymmetric numerical code based on the generalized harmonic formulation 
of the Einstein equations, which is regular at the axis. We test the code by investigating gravitational 
collapse of distributions of complex scalar field in a Kaluza-Klein spacetime. One of the key issues of 
the harmonic formulation is the choice of the gauge source functions, and we conclude that a damped 
wave gauge is remarkably robust in this case. Our preliminary study indicates that evolution of 
regular initial data leads to formation both of black holes with spherical and cylindrical horizon 
topologies. Intriguingly, we find evidence that near threshold for black hole formation the number 
of outcomes proliferates. Specifically, the collapsing matter splits into individual pulses, two of which 
travel in the opposite directions along the compact dimension and one which is ejected radially from 
the axis. Depending on the initial conditions, a curvature singularity develops inside the pulses. 



I. INTRODUCTION 

In general, a detailed investigation of fully nonlinear 
gravitational dynamics is impossible by other than nu- 
merical means. Luckily, the numerical methods have 
recently reached the level of maturity that finally al- 
lows addressing many long-standing puzzles. Perhaps 
the most remarkable is the progress achieved in solving 
a general-relativistic two-body problem — the coalescence 
of black holes [H-Q- Driven by the gravitational wave 
detection prospects, the problem of the collision of two 
black holes or neutron stars continues to be the central 
front where an overwhelming majority of numerical rel- 
ativity research is done. Fortunately, the computational 
methods employed there are portable and — as demon- 
strated below — can readily be applied on other problems 
of interest. 

The success of the numerical simulations was backed 
up by a parallel development of the software and the 
hardware, which provided the necessary computational 
resources. The rapid hardware evolution combined with 
the persisting regularity problems in axial symmetry 
eventually led to a direct transition from highly symmet- 
rical spherical configurations to fully general 3D situa- 
tions without any symmetries at all, essentially bypass- 
ing the intermediate axisymmetric case. However, here 
we argue that important theoretical and practical reasons 
exist to explore axisymmetry better, and we describe a 
new regular numerical code that, we believe, will be ca- 
pable of achieving this. 1 

Before describing any concrete setup we would point 
out one important possible use of an axisymmetric code, 
specifically, that it can be regarded as an efficient "cal- 
ibration tool" for more general 3D codes. 2 Indeed, we 
expect that an intrinsically axisymmetric code applied to, 
say, the head-on collision of two black holes would be ca- 
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1 See [ij- lllll for alternative approaches. 

2 and as a probe of reliability of the cartoon methods 0, Il2ll used 
to effectively simulate axisymmetric spacetimes in 3D. 



pable of following the evolution and the resulting gravita- 
tional radiation more accurately compared to Cartesian 
3D numerical implementations, both because of explicit 
use of the symmetry and since higher numerical resolu- 
tion can be employed for given hardware resources. 

Several interesting open problems arise in axisymmet- 
ric gravitational collapse situations. In particular, it re- 
mains unclear whether or not the weak cosmic censorship 
is violated in collapse of prolate Brill waves [M El- An 
independent observation of universality in critical col- 
lapse of gravitational waves [l4| is pending, as well as 
further investigation of the non-spherical unstable mode 
that apparently shows up at threshold for black hole for- 
mation in axisymmetric collapse of a scalar field [lSj . A 
basic problem of mathematical relativity concerning the 
stability of black holes with respect to nonlinear axisym- 
metric perturbations, can be equivalently addressed in a 
collapse situation by computing how fast a newly formed 
black hole radiates away higher multipole moments. 

In addition, we shall mention other, rarely cited in 
the context of numerical relativity, axisymmetric systems 
which are of great interest in theoretical work on higher- 
dimensional gravity. One of the most basic motivations 
for studying higher-dimensional spacetimes relies on the 
observation that the Einstein equations, describing clas- 
sical General Relativity (GR), are independent of the 
spacetime dimension. Nevertheless, certain properties of 
the solutions to these equations vary dramatically with 
the dimension. A striking example is that axisymmetric 
black holes in dimensions greater than four do not neces- 
sarily have spherical horizons, but also admit horizons of 
toroidal topology [l6j]. Moreover, and in sharp contrast to 
their four-dimensional counterparts, higher dimensional 
black holes do not respect the Kerr limit [ 17[ . and they 
are unstable in certain range of parameters One of 
the fundamental unresolved puzzles [ltl H3] is whether or 
not the instability leads to fragmentation of the horizon 
and exposition of the inner singularity, hence violating 
the cosmic censorship hypothesis. 

Since it is unlikely that any of these these problems 
can be addressed analytically in a systematic manner, 
we turn to computational methods. However, solving 
the Einstein equations numerically is notoriously difficult 
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and depends crucially on the way these equations are 
formulated and evolved. In this paper we focus on the 
generalized harmonic (GH) formulation (2ll - [23l | that has 
recently gained popularity because of its great success in 
the simulations of black hole binaries [ij-yl, 0, H3 ■ 

In a nutshell, the GH approach is a way to write the 
field equations such that the resulting system is mani- 
festly hyperbolic, taking the form of a set of quasilinear 
wave equations for the metric components. As the name 
suggests, the GH method generalizes the harmonic ap- 
proach, which achieves strong hyperbolicity by choosing 
the harmonic spacetime coordinates. In the GH approach 
much of the coordinate freedom is regained through the 
introduction of certain gauge source functions, which also 
maintains the desirable property of strong hyperbolicity 
of the field equations. In fact, the source functions can 
be thought of as representing the coordinate freedom of 
the Einstein equations, and when constructing solutions 
of the equations, via an initial value approach, for exam- 
ple, they must be completely specified in some fashion. 
Choosing the source functions in a controlled way is a key 
issue of the GH technique and after testing several recent 
prescription, we conclude that a damped wave gauge [25[ 
is remarkably robust in collapse situations. 

Applications of the GH approach in spherically- 
symmetric situations were studied in (2(| and here we 
employ the method in axisymmetry. In both cases it is 
natural to use coordinates in which the symmetries of 
the spacetime are explicit. However, these coordinates, 
are formally singular: at the origin in spherical symme- 
try and on the axis in axial symmetry. Thus, the field 
equations have to be regularized in numerical implemen- 
tations; here we describe a regularization procedure that 
is compatible with the GH formulation. 

Our numerical implementation of the GH system is a 
free evolution code that advances initial data by solving 
a set of wave equations. In addition, there are also con- 
straint equations that must be satisfied during the evolu- 
tion. Although the constraints are consistently preserved 
in the GH approach in the continuum limit, in numerical 
computations at finite resolution, constraint violations 
generically develop. In order to maintain stability these 
deviations must be damped and we discuss an effective 
method that achieves this. 

Having in mind implications for higher-dimensional 
GR we test our new code by studying gravitational col- 
lapse of a complex scalar field in a _D-dimensional Kaluza- 
Klein (KK) spacetime. This background, which has a 
single compact extra dimension curled into a circle, is 
a classical example of a higher-dimensional compactificd 
spacetime that in certain limits can appear four dimen- 
sional, for example when the size of the compact dimen- 
sion is small. 3 Assuming spherical symmetry in the infi- 
nite (D — 2)-dimensional portion of the space makes the 



problem 2+1 that depends on time and two spatial co- 
ordinates: one in the radial direction, and one along the 
KK circle. 

We perform a series of numerical simulations where 
the initial distribution of the scalar matter is freely spec- 
ified and the outcome of the evolution depends on the 
"strength" of the initial data as well as on its topology. 
The weak data correspond to the dispersion of relatively 
dilute pulses, while a typical strong data configuration 
leads to black hole formation or nearly does so. Our pre- 
liminary results indicate a wide range of the black hole 
forming scenarios, including how many holes form and 
of what topology. For instance, a static distribution of 
matter centered at the axis and localized in the KK di- 
rection with the energy-density above certain threshold 
collapses to form a black hole with a quasi-cylindrical 
horizon, smeared along the extra-dimension. Data with 
the energy- density below this threshold evolve to form a 
quasi-spherical horizon centered around the initial mat- 
ter distribution. By further reducing the density we find 
that resulting black holes become progressively smaller, 
and at some critical density a pulse of matter is emitted 
radially away from the axis and a curvature singularity 
develops inside it. We find that for slightly lower initial 
densities the evolving matter splits into several pulses 
and two of them individually collapse to form the singu- 
larities moving apart along the KK dimension. Finally, 
when the initial matter distribution becomes dilute below 
a certain limit, no black holes or curvature singularities 
are created. 

In the next section we describe the class of effectively 
2+1 dimensional models where our code can currently be 
applied. In Sec. IIIII we present the basic formulas of the 
GH formulation and discuss the constraints and a method 
to damp their violations. We describe an axis regulariza- 
tion in Sec. IIII Al and boundary conditions in Sec. IIIIBI 
Although, in this work we integrate full _D-dimensional 
equations, in the Appendix we describe an alternative 
approach — one that uses a dimensional reduction on the 
symmetry and integrates the reduced 2+1 equations — 
and compare its performance with ours. Coordinate con- 
ditions and the initial data problem are formulated in 
Sees. IIII CI and IIII Dl respectively. We use several di- 
agnostics to probe the spacetimes that we construct, in- 
cluding computation of asymptotic measurables and ap- 
parent horizons, and describe that in Sec. IIII El After 
elaborating on our numerical algorithm in Sec. IIVI we 
test its performance in Sec. El giving detailed accounts 
of various aspects, such as specific coordinate choices, 
constraint damping, numerical dissipation, and conver- 
gence. We conclude in Sec. IV Dl outline possibilities of 
improvement, and discuss some future prospects. 



II. THE SETUP 



3 In this paper, however, the size of the KK circle is arbitrary. 



We consider a Z?-dimensional spacetime that possesses 
the 0(D — 2) isometry group. We will further assume 
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that the corresponding Killing vectors are orthogonal to 
the closed hypersurface they generate. The symmetry re- 
duces the problem to effectively 2+1, which depends on 
time, t, and two spatial dimensions that we denote by r 
and z. The spatial coordinates can be either infinite or 
finite. For instance, taking D — 4 and assuming asymp- 
totic flatness correspond to the usual four-dimensional 
axially-symmetric situation without angular momentum, 
while setting D — 5 and assuming periodic z and infinite 
r can describe dynamics in a five-dimensional Kaluza- 
Klein background. 

The most general D-dimensional line element with 
these isometries can be written as 



ds = g^v dx^dx v 



g ab dx a dx b + e 2S r 2 dil 2 n . 



(1) 



Here g^ is the D-dimensional metric, dfi 2 is the metric 
on a unit n-sphere, n = D — 3, a, b = 0, 1, 2 running over 
{t, r, z}, and the 3-metric g a b and scalar S are functions 
of t, r and z, alone. 

We take a complex, minimally coupled to gravity scalar 
field to represent the matter of the theory and write the 
total action of the system as 
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V-gD 



Rn — 



dx L 



(2) 



where Gat is the D-dimensional Newton constant. 

In the next section we will describe our strategy to solv- 
ing the equations derived from this action. We will focus 
on asymptotically flat spacetimes times a circle (having 
the topology IR- 0-2 ' 1 x S 1 ) but remark that asymptoti- 
cally de Sitter (dS) and anti-de Sitter (AdS) spacetime 
are also included in our model ([2]) when the potential 
of the scalar field satisfies V(0) — > A with positive and 
negative A, respectively. 



III. GENERALIZED-HARMONIC 
FORMULATION AND ITS REDUCTION TO 2D 
CASE 

In order to numerically solve the Einstein equations 
derived from ([2]), we use the generalized harmonic for- 
mulation. To make the description as self-contained as 
possible, we summarize below basic facts regarding the 
approach, (more details can be found in e.g. @, [H, [26j]) 
and adapt it to the 2+1 situation of interest. 

We begin by noting that whenever isometries are 
present one could perform a Kaluza-Klein reduction on 
them, and in our case ([2]) the reduction yields lower- 
dimensional, 2+1 Einstein equations coupled to the 
scalar, which is related to the size of the n-sphere, and 
the matter. Initially, we had indeed performed such a 



reduction and coded the reduced equations; see the Ap- 
pendix [X] for details and comparison of the methods. 
However, after experimenting with the reduced and the 
full D-dimensional versions of the equations, we found 
that numerical solution of the latter is generically more 
stable. Therefore, in what follows we adopt the unre- 
duced approach. 

The Einstein equations on a Z?-dimensional spacetime 
obtained by varying the action ([2]) can be written in the 
form 



Ruv — SkGnTui/ = 8ttGn T, 



1 



D 



g^T , (3) 



where R^ v is the Ricci tensor, T^ u = d^<& d v )<&* — 
(l/2)g p „(|<9$| 2 + 2 V) is the energy-momentum tensor of 
the matter with trace T. Hereafter we use units in which 
the D-dimensional Newton constant satisfies 8ttGn = 1. 

The Ricci tensor that appears in the left-hand-side 
of ([3]) contains various second derivatives of the met- 
ric components g^ u : these second derivatives collectively 
constitute the principal part of R^ Vl viewed as an op- 
erator on g^y. This principal part can be decomposed 
into a term g a/3 d a pg tiv , plus mixed derivatives of the 
form g ai d atl g lv . Without the mixed derivatives, ^ 
would represent manifestly (and strongly) hyperbolic 
wave equations for the g^ [22j. One can view the GH 
formulation of general relativity as a particular method 
that eliminates the mixed second derivatives appearing 
in ©; see @, i, H El M ■ 

One requires that the coordinates satisfy 



Dx a 



(4) 



where H a = gapH^ are arbitrary "gauge source func- 
tions" which are to be viewed as specified quantities, and 
— r Q = r™ g^ v are the contracted Christoffel symbols. 
One defines the GH constraint 



C a = H a - Ux a . 



(5) 



which clearly must vanish provided (0]) holds, and then 
modifies the Einstein equations as follows: 

R^lv — C(n;v) = T^. (6) 
This last equation can be written more explicitly as 
1 



r^rjL, = 9 (a $ d b) $* + — ^.9a 6 v. 



(7) 



Now, provided that H a are functions of the coordinates 
and the metric only, but not of the metric derivatives — 
namely H a = H a (x >1 ,g) — the field equations ([7]) form 
a manifestly hyperbolic system. The source functions 
H a are arbitrary at this stage and their specification is 
equivalent to choosing the coordinate system ( "fixing the 
gauge"). Determining an effective prescription for the 
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source functions is thus crucial for the efficacy of the 
GH approach, and our strategies for fixing the H a are 
discussed in Sec. IIIICI 

After the coordinates have been chosen, we integrate 
the equations forward in time. Consistency of the scheme 
requires that the GH constraint ([5]) be preserved in time. 
The contracted Bianchi identities guarantee that this is 
indeed the case, since, using those identities, one can 
show 0, [24[ that C a itself satisfies a wave equation, 

ac a + R a u C u = 0. (8) 

Thus, assuming that the evolution is generated from an 
initial hypersurface on which C a — dtC a = 0, and con- 
straint preserving boundary conditions are used during 
the evolution, ([8]) guarantees that C a = for all future 
(or past) times. 

Although the GH constraint is preserved at the contin- 
uum level, in numerical calculations, where equations are 
discrctized on a lattice the constraint cannot be expected 
to hold exactly. It appears that numerical solutions of ([7)1 
can admit "constraint violating modes" , with the result 
that the desired continuum solution is not obtained in the 
limit of vanishing mesh size. However, an effective way of 
preventing the development of such modes in numerical 
calculations exists: one adds terms to the field equations 
that are explicitly designed to damp constraint violations 
(see e.g. [13] )• Following the approach of Pretorius 0, H[ 
that builds on earlier works [28|, |29[ we define the con- 
straint damping terms 

Z^u = K (n^C^ - i.9 p „ n p C^j , (9) 

and solve the modified equations of the form 

Rfiv — C(p;i/) + = T^ v . (10) 

Here, n M is the unit time-like vector normal to the t = 
const, hypersurfaces, that can be written as 

n„ = - (l/<y/=?°) d M t, (11) 

and k is an adjustable parameter that controls the damp- 
ing timescale. Specifically, it is shown in [29| that small 
constraint perturbations about Minkowski background 
decay exponentially with a characteristic timescale of or- 
der n. We note that the constraint damping term con- 
tains only first derivatives of the metric and hence does 
not affect the principal (hyperbolic) part of the equations. 

A. Regularization of the axis, r = 

Having described the GH formulation, we now special- 
ize to the symmetric case. We note first that in our co- 
ordinates ([I} adapted to the symmetry the line element 
of the flat spacetime becomes 

ds 2 = -dt 2 + dr 2 + dz 2 + r 2 dn 2 n , (12) 



and that in this case the source function ((4]) does not 
vanish but becomes 

rrMink T^Mink 
71 

= (0,-,0,(n- 1) cot0i,...,cot0„_i,O), (13) 
r 

where 9i are angular coordinates of the sphere's line el- 
ement dil^. Since near the axis a general spacetime is 
locally flat, the radial component of the source function 
is generically singular at r = 0, diverging as n/r. To 
regularize this radial component, we thus subtract the 
singular background contribution by transforming 

H a ^H a +5lHf mk (14) 

and prescribe gauge conditions using the regular sources. 

Invariance of the line element ([1]) under the reflec- 
tion r — » -r in our case implies that the metric 
components goi and 312 are odd functions of r, while 
goo, 9 11, 522, 502, S and <& are even in r. The GH con- 
straint then dictates that Hi , regularized via (|14p is 
an odd function of r, while Hq and Hi are even in r. 

Moreover, the requirement that the surface area of an 
rt-sphere must vanish at the axis 4 implies gn(t,0,z) — 
exp[2 S(t, 0, z)\. We note that this is an extra condition 
on S, which thus has to satisfy both this relation, as well 
as the constraint that it have vanishing radial derivative 
at r = — specifically that gn— exp(2S') = 0(r 2 ). There- 
fore, at r — we essentially have three conditions on the 
two fields S and gn. In the continuum, and given regular 
initial data, the evolution equations will preserve regular- 
ity, however, in a finite-differencing numerical code this 
will be true only up to discretization errors. As a gen- 
eral rule-of-thumb, the number of boundary conditions 
should be equal to the number of evolved variables in 
order to avoid regularity problems and divergences of a 
numerical implementation. 

An elegant way to deal with this regularity issue in- 
volves definition of a new variable, A, 5 : 

A,^. (15) 

r 

At the axis one then has A ~ 0(r). Therefore, af- 
ter changing variables from 5 to A by using S = 
(1/2) log(<7n — r X) in all equations, and imposing 
X(t, Q,z) — at the axis, one ends up with a system 
where there is no over-constraining due to the demand 
of regularity at r = 0. Crucially, we note that the hyper- 
bolicity of the GH system is not affected by the change 
of variables. 



4 that is, that the radial and areal coordinates coincide at the axis, 
to avoid a conical singularity there. 

5 We note that a similar variable was introduced in [3p| ■ also for 
the purpose of regularization. 
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While we note that a more straightforward regulariza- 
tion method that maintains S as a fundamental dynam- 
ical variable and employs analytical Taylor-series expan- 
sion of the equations in the vicinity of r = can also be 
used in simulations, its reliability degrades in the strong 
field regime, where the regularization that uses A remains 
consistently accurate. 6 



B. The field equations and Kaluza-Klein boundary 
conditions 



With the metric ansatz (fTJ and the regularized source 
function (|14[) . our equations become 8 equations for 8 
variables: six components of the 3- metric g a b, and real 
and complex scalars A and $ correspondingly. Schemat- 
ically the system can be written as 



-9° 9ab,cd — #(a,b) + • 
VA ld rHl, 



D-2 



2 vx+ \0^l. 



D-2 r 
g cd <5>, cd + ... = 0V/0$*. 



9abV, (16) 

(17) 
(18) 



Here ellipses denote terms that may contain the metric 
and it derivatives and/or the source functions, in various 
combinations. These equations are to be evolved forward 
in time starting from the initial (t = 0) time slice, where 
values for the fields and their first time derivatives must 
be prescribed. 

In order to completely specify the problem we have 
to provide boundary conditions which the above equa- 
tions are subject to. In this paper we will be interested 
in a D-dimensional Kaluza-Klein spacetime of topology 
jjD-2,1 x gi^ wnere the z direction is considered peri- 
odic with asymptotic length L namely that z ~ z + L, 
and for the future use we also define half-period L = L/2. 
Asymptotically, this spacetime becomes Minkowski times 
the compact circle and the corresponding boundary con- 
ditions are 



gab -> Vab, A ->■ 0, $ ->■ H a 



Coordinate conditions 



0. 



(19) 



functions H a . The choice that we find to perform best in 
our case is a variant of the damped-wave gauge (DWG) 
condition proposed recently in [25| (see also |32j |) 



F a = 2 pi log ( — ) n a - 2 fx 2 a 1 j a 



(20) 



where y a b = g a b + "a^b is the spatial metric whose de- 
terminant 7 = detjij — (gu g^i — gf 2 ) exp(nS') has the 
factor r™ removed in accordance with the regularization 
(TTJJ, and fii.2 andp are free parameters. 7 Since the gauge 
function (|20[) depends only on the metric, we can simply 
set H a = F a without destroying hyperbolicity of our sys- 
tem. Below we refer to this approach as the algebraic 
DWG condition. 

An alternative method that preserves the hyperbolicity 
was originally devised by Pretorius for binary black hole 
simulations jl|, Q and has also proven useful in the stud- 
ies of the gravitational collapse of scalar field in spherical 
symmetry [26]. This strategy elevates the status of the 
H a to independent dynamical variables that satisfy time- 
dependent partial differential equations. The evolution 
equations for the H a are designed so that the ADM kine- 
matic variables — lapse a and shift f3 l which (implicitly) 
result from the time development — have certain desirable 
properties. For example, the equation for Hq is tailored 
in an attempt to keep the value of the lapse function of 
order unity everywhere — including near the surfaces of 
the black holes — during the evolution. 

One specific prescription for achieving this type of con- 
trol evolves the gauge source functions according to 



a" 



Hi = 0, 



(21) 



where □ is the covariant wave operator, and ao, £x> £2 and 
q are adjustable constants. 8 Thus the temporal source 
function satisfies a wave equation similar to those that 
govern the metric components in the system (1101) . The 
first term on the right-hand-side of (|21l) is designed to 
"drive" H t to a value that results in a lapse that is ap- 
proximately ao- The second, "frictional" term tends to 
confine Ht to this value. For the case of the spatial coor- 
dinates, Pretorius found that the simplest choice of spa- 
tially harmonic gauge, Hi = 0, was sufficient in simu- 
lations of binary black hole collisions. A slight gener- 
alization of this technique was considered in [34| where 
instead of using Hi — 0, the spatial components of the 



As we have already mentioned, fixing the coordinates 
in the GH approach amounts to specifying the source 



See, however, [26l,l3lJ for accurate simulations of the strong grav- 
ity regime in spherical symmetry where Taylor-series approach 
is used. 



We note that the spatial part of DWG is essentially a version of 
the popular T-driver condition |25l |33| . 

In certain situations it is convenient to assume that £1 and £2 
are given functions of space and time rather than mere constants. 
For example, one might require that the gauge driver is switched 
on gradually in time, or that it be active only in certain regions, 
e.g. in the vicinity of a black hole, and that its effect vanish 
asymptotically. 
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source functions are evolved according to 



ft 



(22) 



where £3 is an additional parameter. 

Variants of gau ge drivers were further investigated in 
the recent [25L l26Ll35j] . Specifically in [2^] the following 
hyperbolic first-order drivers were proposed 

d t H a - p%H a = -v (H a - F a ) + W a , 

d t W a + r)W a = -r 1 p i d i H a , (23) 

such that all time-independent solutions of this system 
satisfy H a = F ai where F a are certain predetermined 
target gauge functions, for instance (|2U|) . and the param- 
eters v and r\ are freely specified. 

In Sec. IV Al we compare the performance of these 
strategies and argue that in collapse situations the alge- 
braic DWG approach is the most robust of all. Namely, it 
does not require an extensive fine-tuning of parameters, 
and it is long-term stable. 



D. Initial data 



We will assume initial harmonic coordinates, which im- 
plies, see (dJ), that the lapse can be determined in terms 
of tp as 



a(0, r, z) = ip 2 



(28) 



After substituting this into equation (|2"T)l we solve for ip 
and initialize our basic variables and their derivatives: 



9ij{0,r,z) = ip^Sij, A(0, r, z) = 0, 
d t gij(0,r,z) = d t A(0, r, z) = 0. 



(29) 



We finally note that in order to use the DWG conditions 
correctly, we have to smooth out the transition from ini- 
tially harmonic to later DWG coordinates. As described 
in Sec. IV Al this can be done stably by multiplying F a 
in (|20l) by a time-dependent factor that gradually grows 
from to 1 . For the driver version of DWG we also ini- 
tialize d t W = W = at t = 0. 



E. Spacetime diagnostics 

We employ several diagnostics in order to probe the 
geometry of the spacetimes we construct. 



We now consider specification of initial data, which 
are values for the fields and their first time derivatives 
at t = 0. For simplicity we restrict attention to time- 
symmetric initial conditions for the metric and |$|. 

Given this assumption, initial data for the scalar field 
reduces to the specification of $(0, r, z), which we take 
to have the form of a generalized Gaussian, 



$(0,r,z) = $ e 



-[(l-e 2 r )(r-r a ) 2 + (l-el)(z-z ) 2 }/A 2 



(24) 



where $0 — |^o|exp(i<p) is a complex amplitude and 
ro, Zq, e r , e z , A and cp are real adjustable parameters, sup- 
plemented by the choice 



d t $(0, r,z) = ioj $(0, r, z), 



(25) 



that satisfies 9t|$|t=o = 0, where u> is another parameter. 

The momentum constraint is satisfied for our initial 
data, and writing the initial metric as 

ds 2 = -u 2 dt 2 + i/j 4 (dr 2 + dz 2 + r 2 dnl) , (26) 

one can show that the Hamiltonian constraint becomes 
an elliptic equation for ip(0, r, z) 

a 2 ^ + o 2 ^ + -d r i> + ^1 [{d r 4>) 2 + (a^) 2 ] + 



2(n 



i(l^l 2 + |a z $| 2 + a-Wl<i>l 2 ) + 



0, 



(27) 



that is subject to the boundary conditions, d r ijj\ r —o 
0, d z ip\ z=0 = d z ip\ z = L = and ip\ r ^oo = 1- 



1. Asymptotic charges: Mass and tension. 

Far away from an isolated system a natural radial coor- 
dinate is well defined by comparison with the flat back- 
ground, and the charges that characterize the solution 
can be found from the asymptotic radial behavior of the 
metric functions. In asymptotically flat spacetimes with 
no angular momentum there is exactly one charge: the 
ADM mass. However, in our system with the extra com- 
pact dimension, the tension, associated with varying the 
length of the compact direction, can also be defined. A 
derivation of this result can be found in [36|, but 
here we note why the appearance of an additional charge 
could be expected. Asymptotically the metric becomes 
z-independent since from the lower-dimensional perspec- 
tive the z-dependent Kaluza-Klein modes are massive 
(with the discrete masses m n — 27rn/L, n > 0) and 
decay exponentially ~ cxp(— m n r). In this situation 
asymptotic Kaluza-Klein reduction is possible with the 
result that g zz behaves effectively as a scalar field that 
carries an (unconserved) charge. Designating the asymp- 
totic fall-off of the metric functions gu and g zz as 



f]it 



9z 



= -1 



f 



2a 
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2b 

„D-4 ' 



(30) 



the mass and tension of the solutions are defined as [3t 



rn 
tL 
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(31) 
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FIG. 1: We are locating apparent horizons of two types: one, 
having a spherical topology, is conveniently parametrized by 
p — R(x) (plotted here as centered at (0,0)), and another, 
having a cylindrical topology, is parametrized by r — R(z). 

where 0„ = 27r((™ +1 )/ 2 )/r[(n + l)/2] is the surface area 
of a unit n-sphere, and L is the asymptotic length of the 
KK circle. 



2. Apparent horizons 



Substituting the expressions (|32|) or ([34]) into (|36|) 
yields ordinary second order differential equations for 
R(x) or R(%) correspondingly. The equations could be 
solved by "shooting", subject to appropriate boundary 
conditions, however, here we instead use the following 
point-wise relaxation method, that is more suitable for 
parallel numerical implementations. We start with sup- 
plying an initial guess for the entire function R and iter- 
ate the parabolic equation 

— = -9(R",R',R,g ab ,dg ab ,x a ), (37) 

in unphysical "time" r until a solution is found. This 
equation implies that depending on whether R(t) is in- 
side or outside the horizon at any given moment, it will 
expand or shrink in the next instant. In numerical sim- 
ulations the hope is that the initial R will "flow" to the 
apparent horizon in finite time. We find that in our case 
this indeed happens when the initial guess is "reasonably 
close" to the final solution. 




It is well known that a process that concentrates suffi- 
cient mass-energy within a small enough volume can lead 
to the formation of a black hole. In numerical calcula- 
tions based on a space-plus-time split, black hole forma- 
tion is often inferred by the appearance of apparent hori- 
zons. An apparent horizon is defined as the outermost of 
the marginally trapped surfaces, on which future-directed 
null geodesies have zero divergence. Specifically, in our 
case we will be searching for simply connected horizons 
of two categories, see Fig. [T] (i) Quasi-spherical horizons 
of topology S n , in which case the surface describing the 
horizon can be taken as 



f( P ,x) = P-R(x) 



where 



P = V( r - r o) 2 + {z - z y 
r — 7'o = p sin x 
z- z Q = p cos x, 



(32) 



(33) 



where (ro, Zq) is the point where the horizon is centered, 
and (ii) Horizons of topology x S 1 smeared along 

the z direction, which do not intersect the axis. In this 
case a convenient parametrization is 



f(r,z) = r-R(z). 



(34) 



In either case we define an outward-pointing spacelike 
unit normal to the surface / = 0, 



(35) 



The vanishing of the divergence, 9, of the outgoing null 
rays defined by l a = s a + n a can be expressed as 



h a0 - s a s )VJp = 0. 



(36) 



IV. NUMERICAL APPROACH 

Here we describe our strategy for the numerical solu- 
tion of the GH system (with a scalar matter source) in 
Kaluza-Klein spacetime. 



A. The numerical grid and the algorithm 

We cover the t-r-z space by a discrete lattice denoted 
by (t n ,r iy Zj) — (nAt,iAr,jAz), where n, i and j are 
integers and At, Ar, and Az define the grid spacings 
in the temporal and two spatial directions, respectively. 
As described in the next section, the spatial domain is 
compactified, and hence a grid of finite size N r x N z ex- 
tends from the axis to spatial infinity. Approximations 
to the dynamical fields, collectively denoted here by Y, 
are evaluated at each grid point, yielding the discrete 
unknowns Y^" = Y(t n , rj, zj) = Y(n At, i Ar, j Az). In 
the interior of the domain, the GH equations and the 
gauge-driver equations are discretized using 9 0(h 2 ) finite 
difference approximations (FDAs) , which replace contin- 
uous derivatives with the discrete counterparts, examples 
of which are given in Table H As in 0,0, |2| our scheme 
directly integrates the second-order-in-time equations. 

Following discretization, we thus obtain finite differ- 
ence equations at every mesh point for each dynamical 
variable. Denoting any single such equation as 



= 0. 



(38) 



9 Here h stands for any of the mesh-spacings At, Ar and Az that 
define our numerical grid. 
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dtY 
d r Y 
3 2 t Y 
d 2 r Y 
8lY 
d 2 rz Y 



d r Y 
dlY 



Centered derivatives 



( Y Z +1 ~ K J - 1 )/(2At) 



i 

+1,3 

(\rn -yn 



Oftu - n-i, 3 )/(2Ar) 
(F l "+ 1 -2F l " i + F l 7 1 )/(At) 2 
2Y£+Y£ lii )/(Ar) a 

F^+Y^O/^AiAr) 
y," + i,^i + y"i,. 1 -i)/(4ArA z ) 



V n + 1 



Backward derivatives 



(2 — 5 + 4 y/+ 2,j 



3^+^ 2 ,,)/(2Ar) 



>7k,)/(Arr 



TABLE I: Examples of second order finite-differencing ap- 
proximation to derivatives calculated at a grid point (n,i,j) 
that we use in our discretization scheme in the interior of 
the domain (centered stencil), and at the excision boundary 
(backward stencil) of the numerical grid. 



we then iteratively solve the entire system of algebraic 
equations as follows. 

First, we note that for those variables that are governed 
by equations of motion that are second order in time, our 
0(h 2 ) discretization of the equations of motion results in 
a three level scheme which couples advanced-time un- 
knowns at t n+1 to known values at retarded times t n and 
i n_1 . In order to determine the advanced-time values for 
such variables, we employ a point- wise Newton- Gauss- 
Seidel scheme: starting with a guess for Y™^ 1 (typically, 

we take Y™^~ = Y™A we update the unknown using 



Y, 



yn+l 



Jy 



(39) 



l ,3 



Here 
tion 

m+\ 



TZy is the residual of the finite-difference equa- 
evaluated using the current approximation to 



Y^j , and the diagonal Jacobian element is defined by 



Jy 



hi 



1,3 



dY n+1 ' 



(40) 



In the cases where we used gauge drivers (|2"3")l we found 
that an iteration based on the Crank-Nicholson dis- 
cretization scheme of the corresponding first order equa- 
tions performed well. Specifically, writing any such equa- 
tion schematically as Y = fy(Y, BY, . . .), we update using 



Y 



• ■j 



y» 



1 



n+1 



+ h 



(41) 



We iterate (|3"9"f and (dJ) over all equations until the total 
residual norm, see (|49|) . falls below a desired threshold. 

In order to inhibit high-frequency 10 instabilities which 
often plague FDA equations, our scheme incorporates 



10 "High-frequency" refers to modes having a wavelength of order 
of the mesh spacings, Ar and Az. 



explicit numerical dissipation of the Kreiss-Oliger (KO) 
type [3a |. Following 0|, at the interior grid points, 
{(i,j)\ 2 < i < N r - 2, 2 < j < N z - 2}, and for 
each dynamical variable we apply a low-pass KO filter 
by making the replacement 



Y 



dij — 



Y, 



26 (^ _2 'J — ^Y-i,j — 41i+i,j + Yi + 2j + 
+ Yi t j-2 — ^Yi.j-i — 4Yij + i + Yij+2 + 



12 Y 



(42) 



at both the t™ -1 and t n time-levels before updating the 
t n+l unknowns. Here cko is a positive parameter satisfy- 
ing < €ko < 1 that controls the amount of dissipation. 
An extension of the dissipation to the boundaries [2| was 
also tried, but this has not resulted in any significant 
improvement of the performance of the code. 



B. Coordinates and boundary conditions 

While the physical, asymptotically flat (times a circle, 
in our case) spacetime extends to spatial infinity, in a nu- 
merical code one can only use grids of finite size. Instead 
of using a standard strategy that deals with this issue 
by introducing an outer boundary at some finite radius 
where approximate boundary conditions are imposed, we 
adopt another technique — proven successful in previous 
work in numerical relativity, see e.g. 0,[2(J[26[ — and com- 
pactify the spatial domain. We find that compactifying 
the radial direction and imposing the (exact) Dirichlct 
conditions (|T9"|) at the edge of the domain works well, 
provided that we use sufficient dissipation. In particu- 
lar, it is known that due to the loss of resolution near 
the compactified outer boundary (assuming a fixed mesh 
spacing in the compactified coordinates) , outgoing waves 
generated by the dynamics in the interior will be partially 
reflected as they propagate toward the edge of the com- 
putational domain, and these reflections will then tend 
to corrupt the interior solution. By adding sufficient dis- 
sipation one can damp the waves in the outer region, 
attenuating any unphysical influx of radiation. This en- 
ables one to use the compactification meaningfully. 

The results presented in this paper were obtained using 
the compactification of the form 



1 



(43) 



where the compactified f ranges from to 1 for values 
of the original radial coordinate r £ [0, oo). In practice, 
we define a uniform grid in the compactified f and use 
chain-rule to replace derivatives with respect to r with 
the derivatives with respect to f in all dynamical equa- 
tions. The asymptotic boundary conditions (|19[) at f — 1 
are then imposed exactly: g a b = rj ab ' , A = 0, $ = 0, and 
H a = W a = 0. 
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We have previously described the boundary (regular- 
ity) conditions at f = r = in Sec. IIII Al Denoting by 



Y, 



n+l 



x j , for j = 2, . . . , N z — 1, the advanced-time value at 
the axis for any of the variables, goo, gu, g22, go2,Ho, and 
Hi that have vanishing derivative at r — 0, we use the 
update Y]™/ 1 = (4F™+ 1 - Y$ l )/Z, which is based on 
a second-order backwards difference approximation (see 
Table of d r Y — d?Y — 0. For the quantities 501, 312, A 
and H\, which are odd in r as r — ¥ 0, we simply use 

Kt = °- 

For simplicity, in this paper we consider only configu- 
rations that have reflection symmetry about z = which 
together with periodicity in z, implies 



d z Y\ z=0 = d z Y\ z=L = 0. 



(44) 



We update points at z — and z — L using backward 
difference approximation similar to that in Tab. |H In 
particular for i = 2, . . . , N r - 1, we use = (4 F/ l 2 +1 - 

Y^)/3 at z = 0, and = (4^_i - *$"i 2 )/3 

at z = i. 



C. The elliptic equation of initial data 

It turns out that generally the equation (j2"T|) for the 
conformal factor is ill-posed, namely that the linear equa- 
tion governing small perturbations about a solution of 
(|27| does not admit a unique solution for given bound- 
ary conditions 39]. Therefore, an attempt to solve ([27)1 
using standard relaxation methods will generically fail, 
as the relaxation is not guaranteed to converge. How- 
ever, a method circumventing this difficulty exists [39l |. 
and this is through a rescaling $ = $ V> s that transforms 
the Hamiltonian constraint ([27)1 into 



6^ + <9 Z 2 V - 

s 2^,2( s -l/2) 



-d r -0 H 

r ip 



1 



[(9 r V) 2 + (c^) 2 ] + 



4(n + l) 

sip 2 



4(n 

^2(s+l/2) 



i$i 2 [(9 r ^) 2 + (9 z v) 2 ; 

- (9 r |$| 2 a^ + 9 2 |d| 2 9 z v 



+ 



(45) 



$1 



4 (n + l) 

w 2_02( S -2n+5/2) 
4(n + 1) 



l$l 



V> 5 



2 (n + l) 



V s ) = 0. 



By choosing the power s such that terms in this equation 
which are proportional to tl> Pi have only non positive pi 's 
one renders the problem well posed 39]. In the case of 
a free scalar field having the potential V cx I*! 2 , which 
we assume here, this implies that s < —5/2, see a related 
discussion in |40| . 

Note that the physical data are $, not <t, and therefore 
we solve equation (|45p in an iterative manner where we 
start with an initial guess ip = ipo(r, z) and at each itera- 
tion i > 1 update = $ i/j^ s that is then used in P5)l 



in order to solve for ipi+x- Most of the results obtained in 
this paper were generated using ipo — 1 and s = —3. In- 
terestingly, it turns out that provided the initial guess is 
not too distant from the solution, the original equation 
([27]) is numerically stable without the rescaling. This 
happens, for instance, in the weak field regime where the 
guess ^0 = 1 relaxes to the solution without any trouble. 



D. Excision 

We use excision to dynamically exclude from the com- 
putational domain a region interior to the apparent hori- 
zon that would eventually contain the black hole singular- 
ity. This approach relies on the observation that in space- 
times that satisfy the null energy condition and assuming 
the cosmic censorship holds, the apparent horizon is con- 
tained within the event horizon. This ensures that the 
excluded region is causally disconnected from the rest of 
the domain (see [4l| and the references therein for further 
discussion). Operationally, once an apparent horizon is 
found, we introduce an excision surface, Rex, contained 
within the apparent horizon, and such that all charac- 
teristics at R = Rex are pointing inwards. This specific 
characteristic structure eliminates the need for boundary 
conditions at Rex- rather, advanced-time unknowns lo- 
cated on the excision surface are computed using finite 
difference approximations to the interior evolution equa- 
tions, but where centered difference formulas are replaced 
with the appropriate one-sided expressions given in Table 
III 

Currently our apparent horizon finder is only capable 
of locating horizons of the shapes depicted in Fig. [T] 
In this case the radius of the excision surface typically 
satisfies Rex ~ 0.7i?AH, where R is the coordinate radius. 



V. PERFORMANCE OF THE CODE 

In this section we investigate the performance of the 
code in series of simulations that evolve regular ini- 
tial distribution of complex scalar field in five and six 
dimensional Kaluza-Klein spacetime. The code uses 
pamr/amrd infrastructure [42( where our suitably inter- 
faced numerical routines are called by the amrd driver. 
All our results are generated using an initial scalar field 
profile of the Gaussian form (f2"4"|). with fixed values 
ro = zo — and A = 0.25, so that the scalar pulse is 
always initially centered at the axis, see Fig. [21 In ad- 
dition, we set e r = and use e z = 0, unless otherwise 
specified. We choose the initial frequency in ([2"5]) to be 
u) = 20, and the asymptotic size of the KK circle, L = 2. 
The scalar field potential used here is taken to be of the 
form V(§) = m| l^l 2 , where without much loss of gen- 
erality we set m$ = 1. 

Because we mostly use a time-explicit finite difference 
scheme, we expect restrictions on the ratio Ac = Ai/h 
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FIG. 2: We choose the initial distribution of the scalar field 
as a generalized Gaussian (|24[) and assume that the initial 
time derivative is given by (|25p . We mostly consider localized 
initial pulses (such as the one shown here), that are obtained 
by setting e r — e z = in (|24[) . The radial direction is com- 
pactified, and the hypersurfaces z = 1 and z = — 1 identified 
in a Kaluza-Klein spacetime under consideration. 

(the Courant factor) that can be used while maintain- 
ing numerical stability. In the results discussed below 
we chose 0.2 < Ac < 0.5 in the weak-field regime, and 
0.05 < Ac < 0.25 for evolving the strong data. Taking 
larger Ac usually leads to amplification and dominance 
of numerical errors near the axis, that shows up as di- 
verging high frequency oscillations. In most cases we use 
uniform grids of the same size in both spatial directions, 
h = Ar = Az, and our lowest and highest resolution 
simulations have h = 1/16 and h = 1/256, respectively. 
The highest-resolution runs generally required Ac < 0.2 
for stability. 

As expected, the outcome of the collapse depends on 
the strength of the initial data. Low density distri- 
butions describe weakly gravitating scalar pulses which 
completely disperse in all instances. Strong data gener- 
ate spacetimes in which black holes form, or almost form. 
For fixed A, e r , e z ,m^, L and ui a single free parameter 
that controls the strength of the pulse — and hence, the 
outcome of the evolution — is the initial amplitude $o- 
In this case "strong initial data" will refer to situations 
when an increase of the initial amplitude by less than 
10% leads to formation of a curvature singularity, and 
the other data will be called "weak" . The critical ampli- 
tude is at the threshold for the singularity formation. 

For each spacetime that we construct, the asymp- 
totic mass and tension are computed using (I31|) . where 
the constants a and b are found by fitting the metric 
components gu and g zz with the functions of the form 
<?( r ) — <?oo + 9i/r n + gijr n+x in the asymptotic region. 
The errors in the constants are determined by the fitting 
uncertainty which in our case is ~ 1 — 2% (larger for 
weaker data). We find that our initial data sets usually 
have a ~ ni, and it follows from (I31[) that r ~ 0, within 
the numerical accuracy of our method. 



t=0 t ~3 




FIG. 3: Snapshots of the evolution of energy-momentum den- 
sity p = n a n b T a b in the simulation of the initial data defined 
by $o = 0.35 (1 — i) in 5D, using algebraic DWG. The hor- 
izontal axis represents the compactified radial direction, and 
the vertical axis is along the periodic z-direction. The lapse 
of time is measured in units of the total mass. The initial 
collapse of the pulse around the center (r, z) = (0, 0) is ac- 
companied by a growth of the amplitude of p by about an 
order of magnitude. The pulse then bounces off and sev- 
eral outgoing shells form. The shells propagate along the 
periodic z-direction and interact, creating a typical interfer- 
ence pattern. As expected, the z-dependence vanishes in the 
asymptotic regions. In the later stages of the evolution p is 
small, and the frequency spectrum of the oscillations along z 
is dominated by a few lowest eigenfrequencies associated with 
the compact KK circle. 



A useful quantity that illustrates energy-momentum 
distribution is the density p — T a i,n a n b . Figure [3] de- 
picts p at several moments during the 5D evolution of 
weak initial data defined by <&o = 0.35 (1 — i), using al- 
gebraic damped wave gauge (DWG), see (|2"0"jl . Figure 
[3] shows that in the early stages of the evolution the 
pulse (initially localized around the center (r, z) = (0, 0)) 
undergoes a gravitational collapse that is accompanied 
by a growth of the central energy-density. At a later 
time, however, the pulse bounces off while forming sev- 
eral shells, and disperses. The distribution of the energy 
density is anisotropic inside the shells, with a higher con- 
centration of energy occurring along the axis, r = 0, and 
at the equator, z = 0. The matter waves that travel 
along the compact KK circle collide to form a typical in- 
terference picture (see Fig. [3]) and the pattern becomes 
increasingly complicated in the course of time, when more 
and more waves undergo interaction. As expected in the 
KK background the z-dependence of solutions vanishes 
in the asymptotic region, since the z-dependent modes 
are massive and fall off exponentially fast. 
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FIG. 4: The logarithm of the central energy-density, p, in a 
five-dimensional weak field simulation defined by |$o| = 0.55 
that uses algebraic DWG. The resulting spacetime has the 
mass M = 0.67 ± 0.03 and the tension is zero within the 
numerical accuracy. In the early stages of the evolution, the 
pulse collapses and p grows. Subsequently, the energy density 
decreases following a power-law dependence during the period 
lasting from t ~ O(10)M until t ~ O(100)M. After that p 
decays exponentially exp[— i/(800Af)]. The high-frequency 
oscillations are associated with the eigenmodes of the compact 
KK circle. 
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FIG. 5: A log-log plot of the Ricci scalar (in units of M 2 ) as 
a function of the proper time, t pr o P = J ' a dt, both evaluated 
at the center, in several simulations in 5D. The first peak in 
the curves corresponds to the initial collapse phase, and for 
stronger data the Ricci scalar diverges in finite time, signaling 
a curvature singularity (case |$o| = 1.2021). In subcritical 
cases, the curvature decreases after the initial growth, and the 
secondary peaks correspond to collisions of the shells traveling 
along the KK circle and arriving back at (0, 0). 



In order to obtain more quantitative insight into the 
process, we plot in Fig. [4] the evolution of the loga- 
rithm of p at the center. During the highly dynamical 
early epoch, lasting until ~ O(10)M, the field collapses, 
bounces off the center, spreads along the z-direction, and 
starts dispersing to infinity. In the first stage of the dis- 
persion, lasting until ~ O(100)M, the decay of p follows 
a power-law, and in the late times the decay is exponen- 
tial with a characteristic time scale of a few hundreds 
of M . The maximal scalar curvature is attained during 
the initial collapse phase; in the shown simulation the 
curvature reaches the magnitude of order of ~ 0(100) 
(in units of M~ 2 ). The discrete spectrum of the high- 
frequency oscillations consists of the normal frequencies 
fn ~ Ti/L, n = 1, 2, . . . defined by the size of the KK cir- 
cle. We observe that the higher frequency modes decay 
faster, and the late-time spectrum is dominated by a few 
lowest frequency modes, see also bottom right panel of 
Fig. El 

When the initial amplitude of the scalar field in- 
creases, the maximal curvature achieved during the col- 
lapse grows, and when |<j> | surpasses a certain threshold, 
the curvature diverges, signaling the appearance of a sin- 
gularity; see Fig. [5] Currently we are able to estimate 
the critical amplitude with a relatively low precession 
of approximately 1 part in 800. In five dimensions the 
amplitude is between 1.0946 and 1.096 such that initial 
data with |$ | < 1.0946 completely disperses, and the 
data |$o| > 1.096 gives rise to curvature singularities 
and black holes. 

A covariant way to illustrate the distribution of mat- 
ter and the geometry of the spacetime is provided by 



the Ricci 11 scalar, and Fig. [6] shows its evolution com- 
puted in the evolution of our strongest subcritical data 
set, |$ | = |0.774(1 - i)\ ~ 1.0946. The initial pulse col- 
lapses, bounces off and forms several shells, that start dis- 
persing after t ~ 0.2M . While the anisotropy of the en- 
ergy distribution inside the shells is small for weak data, 
it is obvious in the strong field regime. Some matter 
is ejected along the equator, z = 0, and two distinctive 
pulses, moving in the opposite directions, form at the 
axis. The pulses collide and interact, which results in 
smearing of the energy-density along the axis, see bottom 
right panel of Fig. [SI The distribution is not stationary, 
rather the matter is continuously leaking to infinity, and 
the late-time behavior is qualitatively similar to the weak 
field case, shown in Fig. [3] 

It turns out that for higher initial amplitudes, 
| <3>o I ~ 1-096, the pulses that form at the axis are able 
to individually collapse and develop curvature singulari- 
ties, indicated by diverging scalar curvature and energy- 
momentum density while the lapse remains finite (of or- 
der one) everywhere. Presently, our horizon-finder is not 
designed to locate the moving apparent horizons that 
may arise in this case around the singularities, and it 
would be very interesting to verify whether or not such 
engulfing horizons indeed form. The time when the cur- 
vature singularities appear depends on the strength of the 
initial data in the manner that the closer to threshold we 
are the later into the evolution the curvature diverges. 



In this and other figures Ricci scalar is measured in units of M 2 . 
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FIG. 6: The Ricci scalar at several instants of the evolution of 
the initial data characterized by $0 = 0.774(1 — i), currently 
our strongest subcritical data set. The numerical resolution 
used in this simulation is 129 x 129. After the initial collapse 
stage an expanding shell of matter forms. The energy-density 
distribution inside the shell quickly becomes anisotropic with 
most of the energy localized within individual pulses traveling 
along the axis and inside the pulse emitted radially along the 
equator, z — 0. The pulses at the axis undergo interactions 
and spread the energy-density along r = 0. In late stages of 
the evolution, all the matter is radiated away to infinity. 



For instance, in the case of |$o| = 1.096 the pulses, col- 
lapse and become singular as they cross the circle and 
collide at z ~ 1, see Fig [7] However, for |$ | = 1-0975, 
the curvature inside each pulse blows up earlier, when 
they reach z ~ 1/2. 

Intriguingly, we observe that for even stronger initial 
data formation of the curvature singularity ensues differ- 
ently Specifically, the pulse of matter which is usually 
emitted during the initial collapse-bounce stage outwards 
along the equator, is now seen to also be able to collapse 
and develop a singularity. Snapshots of the process are 
shown in Fig. [5]that was obtained in the evolution of the 
data defined by |$o| = 1-099. Since all our attempts to 
detect an apparent horizon, engulfing the singular region 
and the center, have failed we believe that the horizon, if 
it forms in this case, must be localized around the moving 
curvature singularity. 

In Fig. IHlwe plot the Ricci scalar along the equator at 
t ~ 0.37M, shortly before the appearance of a singularity 
at f ~ 0.21 causes the code to break down. Figure 
shows that the maximal value of the Ricci scalar along 
this slice is determined by the numerical resolution: the 



FIG. 7: Snapshots of the Ricci scalar at several instants of the 
evolution of the initial data characterized by $o = 0.775(1 — i) 
at the resolution of 257 x 257. The matter distributes 
anisoropically inside the bouncing shells such that two local- 
ized pulses, traveling along the axis, form. When the pulses 
collide at z = 1 a curvature singularity appears. 



finer meshes we use, the larger the values attained. This 
signals emergence of a curvature rather than a coordinate 
singularity. 

Figure [TU] shows snapshots of the energy- density dis- 
tribution obtained in a supercritical simulation initiated 
by $o = 0.8(1 — i). The resulting spacetime has the total 
mass M = 2 AO ± 0.06 and the tension r ~ 0. The initial 
stages of this evolution are qualitatively similar to those 
of weaker data, however, after the first bounce off the 
center (at t ~ 0.19M) the pulse recollapses and a black 
hole forms as indicated by the appearance of an apparent 
horizon. The energy-density p and scalar curvature both 
blow up in the vicinity of (r, z) — (0, 0) in finite time, 
while the lapse and shift remain finite. 

In order to locate apparent horizons we solve (1371) it- 
eratively starting with some initial guess for the entire 
function that parametrizes the horizon. In the case of 
$ = 0.8(1 — i) we are searching for a horizon that has 
spherical topology. We use the parametrization (|32|) with 
ro = zq = and solve (|57)l initialized by R = 0.1. Set- 
ting the target accuracy to ~ 10~ 3 and using 100 grid 
points to represent the horizon, we are able to solve the 
equation in ~ 1000 iterations. The resulting horizon is 
shown in bottom panels in Fig. [TTJ1 

After an apparent horizon is found we use excision to 
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FIG. 8: The Ricci scalar at several instants during the evo- 
lution defined by $0 = 0.777(1 — i). A curvature singularity 
forms along the equator, inside the lump of matter emitted in 
the radial direction during the initial collapse- bounce process. 
This simulation uses the resolution of 257 x 257. 
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FIG. 9: The Ricci scalar (in units of M~ 2 ) along the equa- 
tor at t ~ 0.37Af, shortly before the evolution defined by 
$0 = 0.777(1 — i) develops a singularity at r ~ 0.21. The 
maximal value that the Ricci scalar attains in this evolution 
depends on the resolution (here h — 1/256), being larger on 
finer numerical meshes, which indicates the geometrical na- 
ture of the emerging singularity. 



remove a region inside the horizon that contains curva- 
ture singularity. However, it turns out that regardless 



FIG. 10: The distribution of energy- density in the evolution 
of the initial data, defined by $0 = 0.8 (1 — i) in a simulation 
that uses the resolution of 201 x 201 . A strong gravitational 
interaction in this case precludes any significant amount of 
matter from escaping to infinity and leads to formation of 
a black hole soon after the first bounce. The black hole is 
centered around (0, 0) and is signaled by the appearance of 
an apparent horizon (designated by a thick dashed line) and 
an explosive growth of the scalar curvature and the energy- 
density. 



of our specific gauge choice the lapse keeps evolving in- 
side the horizon, and continuously decreases until it is 
reaching the magnitudes of order ~ 10 -5 near the ex- 
cision boundary. In this situation truncation errors in 
quantities near Rex occasionally cause the computation 
of non-positive values for the lapse, which immediately 
leads to code failure. Unfortunately, this happens when 
the apparent horizon is still fairly dynamical and con- 
tinues to change its shape and size. Therefore, we are 
currently not able to determine the stationary black hole 
state in the end of the evolution. 

Having described the low mass configurations we will 
now briefly discuss more massive solutions of certain 
type. One initial configuration that we have evolved 
consisted of a nearly uniform along z distribution of the 
matter — achieved by setting e z — 0.995 in (fM|) — with the 
initial amplitude of $0 = 1-2(1 — i). The Ricci scalar and 
the energy-density computed in this evolution are shown 
in Fig. [Tl] together with the resulting apparent horizon 
that appears at t ~ 0.55 and has a cylindrical topology. 
In this simulation we used the algebraic DWG condition, 
and in this case we were not able to locate a surface on 
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FIG. 11: Distributions of the energy-density (on the right) 
and the Ricci scalar (on the left) computed in the evolution of 
a nearly uniform along z initial data defined by $o = 1.2 (1 — 
i), and e r = 0, e z = 0.995, that uses the resolution of 101 x 33. 
A strong gravitational interaction leads to a rapid formation 
of a black string, signaled by the appearance of an apparent 
horizon of cylindrical topology (shown as a thick dashed line). 
Note, that the matter and the curvature inside the horizon are 
strongly localized around the center r — z = 0. 



which all characteristics will point inwards. Therefore, no 
excision was employed in this evolution, and the eventual 
failure of the code was caused by collapse of the lapse at 
the axis. 

While the total mass of this spacetime, M ~ 24.5, is 
determined by the asymptotic fall-off (|3ip. the mass of 
the "black string" (a black hole with the smeared hori- 
zon) can be estimated from the average size of its horizon. 
The last moment of the evolution before the code had 
crashed is shown in bottom panels of Fig. 1111 The hori- 
zon is nearly uniform along z and is located at ?ah ~ 0.3, 
which yields the mass of Mbs = 0-5Rah/(Gn/L) ~ 12, 
where Rah ~ 0.49 is the uncompactified average areal 
radius of the horizon. This value is below the critical 
mass, M c ~ 14, needed for stability of the extended solu- 
tions of this type 18[ . Therefore, at this stage the system 
is probably far from approaching stationarity. Since the 
total available mass is higher than M c several scenarios 
for subsequent dynamics can be imagined. If the black 
hole accretes enough matter and increases its mass above 
critical, it can reach the uniform black-string end-state. 
Otherwise, the evolution would probably have to proceed 
via forming a localized black hole in the first place. This 



black hole then may or may not accrete additional mat- 
ter, and as a result either to grow and become a black 
string or to settle down to a stationary solution of spher- 
ical topology. While further investigation of the process 
is clearly needed, Fig. [TT] shows a developing progres- 
sive localization of energy- density and curvature around 
(0, 0), indicating that formation of a localized black hole 
first is more probable in this specific case. 

We turn now to a detailed description of the perfor- 
mance of the code. Since we have implemented a free 
evolution scheme, we can assess the convergence of our 
numerical solutions by monitoring discrete versions of 
the Hamiltonian and momentum constraints, which are 
defined by contracting the Einstein equations with the 
unit normal vector to the t = const hypersurfaces, i.e. 
M a = n a (G a b — T a i,), where G a b is the Einstein tensor. 
One way of doing this involves evaluation of the following 
L2-norm of finite-differencing variables 



\Y\\l 2 = 



1 



N r N, 



N r ,N z 

E 



1/2 



(46) 



In the next section, we compute the L2-norms of the 
constraints at each time step and examine their behavior 
as a function of various parameters of the problem. 

A. Coordinate conditions 

In the weak-gravity regime where an initial pulse of 
matter completely disperses to infinity we find clear ad- 
vantage of the algebraic DWG conditions. They are 
robust, almost do not require fine-tuning, and are ex- 
tremely stable. Additionally, as described in the next 
section, the constraints remain well preserved in this case, 
such that only very small damping (controlled by the pa- 
rameter k, see (O) needs to be added to the equations. 

Since at t = 0, we assume harmonic conditions H a = 
0, choosing the source functions according to (|2"0|) at 
t > will create discontinuity in the temporal compo- 
nent of the gauge source function. Therefore, we multi- 
ply the sources F a in (|2T))) by a time-dependent function 
(1 — exp(— i/to))(l + s exp(— t/ttj), where to,ti, and s 
are parameters. While the first factor is used to gradu- 
ally turn the sources on, the second factor is introduced 
in order to improve late time stability. For instance, 
while we found that weak data defined by |$o| ^5 0.35 
can be simulated using to ~ — 1 and s — 0, the pa- 
rameters to ~ — 0.2, s ~ 0.1 — 0.5 and t\ ~ 10 were 
required in order to maintain regularity during evolution 
of stronger data. In addition, we found that long-term 
stability of the strong data evolutions improves when the 
factor exp(n S) is removed from the definition of the de- 
terminant 7 = (gn gi2 — ffL) 1 ^ 2 exp(n S) used in (j2"U)) . 

The simulations that employ algebraic DWG are sta- 
ble for a range of the parameters [i\ and [i2- Basically, 
any values of these parameters of order one can be suc- 
cessfully used in the collapse situations. Nevertheless, 
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FIG. 12: While evolutions are usually stable for any values of 
order one of the parameters /ii,2 in (|20[) . instabilities develop 
if these parameters are too small or too large. Here we show 
the behavior of the norms of the constraints (|46[) in the simu- 
lation of the date defined by $o = 0.3 in 6D. The evolution is 
stable for 0.3 \i\ t i i$ 5 and the decay rate of the constraint 
violations is optimized for |tii,2 — 2. 



we still find that certain choices of fi\.2 perform better 
than others and preserve the constraints with a greater 
accuracy. This is illustrated in Fig. [12] where we show 
norms (|46|) of the constraints as a function of time. Al- 
though, the norms decrease in all the instances, the con- 
straints violations are the smallest for = //2 — 2. 

After experimenting with the driver version of DWG 
condition (|23|) . we find that it performs comparably to 
the algebraic DWG in early stages of the evolution until 
t ~ 10M. However, late-time behavior is fairly sensitive 
to the driver parameters and generically develops coor- 
dinate singularities on a timescale t ~ 30 — 100M. A 
considerable constraint damping was usually required in 
these cases. In comparison, a DWG-driver simulation 
with parameters identical to the algebraic DWG evolu- 
tion shown in Fig. @] had required at least 50 times 
stronger damping in order not to diverge immediately, 
but even then the evolution remained stable only until 
t ~ 100M. 

The performance of the gauge drivers (|2"T1) and (|22l) is 
similar to that. It was already reported in 26] that those 
drivers — known to perform well in Cartesian implemen- 
tations [J, H| — are considerably less stable in spherical 
symmetry, and here we notice the same. Specifically, we 
find that pure harmonic coordinates are useful only for 
simulating the dynamics of very weak initial data with 
the scalar amplitude |$o| ~ 10 -3 (and the corresponding 
mass M IS 10~ 4 ). For larger values of |$o| the lapse func- 
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FIG. 13: The behavior of the Hamiltonian and the momentum 
constraints, Al a = n a (G a b — T a b), as a function of the dissipa- 
tion parameter tKO in a 6D simulation with $o = 0.3 (1 — i) 
and the resolution 33 x 33. While for €ko ~ 0.04 the con- 
straints blow up and the code diverges, above this threshold 
the evolution is stable for arbitrary long periods of time, and 
any constraint violations decrease in the course of evolution. 

tion collapses in the locus of maximal matter concentra- 
tion and coordinate singularity forms. Although we were 
able to delay the pathology by using of one the drivers 
(|21l22p , in no instance was it possible to eliminate it com- 
pletely. Generically, the coordinates in these simulations 
become singular after a time of approximately a few tens 
of M. An extensive fine-tuning of the driver parameters 
together with stronger constraint damping and the nu- 
merical dissipation enables one to extend somewhat the 
duration of the regular evolution, up to t ~ 100M. For 
example, a parameter setting that kept the evolution of 
the | <E>o I = 0.5 data regular until t ~ 100M consisted of 
e KO = 0.34, A c = 0.1, & = 0.6, & = 0.8, & = 0, a = 
1, q = 3 and N r x N z = 41 x 41. 

Our preliminary experiments in supercritical regimes 
indicate that the algebraic DWG still outperforms the 
driver conditions. While in all our simulations where a 
black hole forms the coordinates eventually become sin- 
gular near the excision surface, the simulations employing 
algebraic DWG remain stable for longer. It is presently 
unclear whether the coordinate pathology is caused by 
the strong gravitational field or by the introduction of 
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FIG. 14: The lapse at the center and the energy-density 
averaged along the axis are shown as functions of time for 
various values of eno in a subcritical 5D simulation, with 
4>o = 0.4(1 — i). Provided we use sufficient dissipation, 
(zko > 0.08 in this case), the behavior of dynamical vari- 
ables is essentially independent of the specific exo in early 
time until ~ 900 M, see the upper inset in the top panel. 
While at the late times the absolute values of variables in the 
simulations with stronger dissipation are generically smaller 
than those in simulations using less dissipation, the variations 
of the fields are very small at this stage in all instances, see 
bottom inset. 



the excision surface. More experiments are required and 
will be reported elsewhere. 



B. Dissipation and constraint damping 

The explicit numerical dissipation that we add to our 
scheme in (|42[) is an important ingredient affecting the 
long-term stability, and Fig. fT~3l illustrates this. There we 
fix $ = 0.3 (1 - i) exp[-(r 2 + z 2 )/0.25 2 ], use resolution 
of 33 x 33 and algebraic DWG conditions with /zi,2 = 2 
and q = 1/2 in a 6D simulation. We find that without 
the dissipation the constraints quickly diverge. However, 
for tKO above a certain minimal value (cko > 0.04 in 
this case) the evolution stabilizes. The specific threshold 
value increases slightly when denser grids and stronger 
initial data are used; however, taking €ko — 0.15 — 0.4 
was usually a safe choice for the simulations described in 
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FIG. 15: The dynamics of the norm (T4U)) of the GH con- 
straints (0) depends on the amount of the damping. In sim- 
ulations of the initial data defined by $(0,r, z) — 0.3(1 — 
i) exp[— (r 2 + z 2 )/0.25 2 ] that employ algebraic DWG any con- 
straint violations decrease in time, provided k < 0.0185A -1 ; 
for larger values of k the constraint violations diverge. An op- 
timal preservation of the constraints in this case is achieved 
for k ~ 0.05/(1 + 0.05 1). 



this paper. 

In Fig. [14] we depict the lapse a(t, 0, 0) and the mat- 
ter energy-density p(t, 0, z), averaged along the axis, as 
functions of time for several choices of £ko m 5D sim- 
ulations of weak scalar pulse with $o = 0.4 (1 — i) us- 
ing algebraic DWG. The simulations converge for a wide 
range of the values of the dissipation parameter, provided 
£a'o ~ 0.08. Figure [T4l demonstrates that in this case the 
specific values of cko have only a marginal effect on the 
early dynamics, until approximately 1000 M. However, 
after that time the variables computed in the simulations 
using unequal dissipation parameters begin to differ, and 
stronger dissipation generically implies smaller late-time 
amplitudes. The absolute values of the amplitudes are 
usually small at this stage (typically below ~ 10~ 4 ). 

Another factor influencing stability is the constraint 
damping term ([5]), which we add to the Einstein equa- 
tions in (fTUj) . We find that quite generically the long-term 
stability improves when the damping of the constraints 
vanishes in the asymptotic regions. For this reason we 
multiply k by a factor (Rq/R) in the regions where the 
areal radius satisfies R > Ro for some large Ro (typically 
Rq ~ 20) in order to gradually turn the damping off. 

Figures [T5] and [T|5] illustrate the effect of the damp- 
ing in the case $(0,r, z) = 0.3(1 — i)exp[— (r 2 + 
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FIG. 16: Same simulation and same notations as in Fig |15l 
The Hamiltonian and the momentum constraints are asymp- 
totically satisfied provided a correct amount of damping is 
used. 



z 2 )/0.25 2 ], cko — 0.5 in 6D simulations that use alge- 
braic DWG with ni t 2 = 1 and q = 1/2, and the resolution 
of 33 x 33. In this configuration the evolution is stable for 
rather small values of k including that of K — 0; however, 
the quickest asymptotic decrease of the constraint viola- 
tions is achieved for k = 0.05/(1 + 0.05 1). For the values 
of k greater than a certain value — k > 0.25 in this spe- 
cific simulation — the code diverges. The time-dependent 
factor 1/(1 + 0.05 1) is less crucial in shorter subcritical 
simulations, but it helps to improve stability on the long 
time scales of t > 2000 M. 

Interestingly, typical values of the damping parame- 
ter k in algebraic DWG evolution are small compared 
to the inverse of any typical length scales of the prob- 
lem (set e.g. by the initial width of the scalar pulse, A, 
by the size of the KK circle or by the size of the black 
hole). Moreover, certain weaker initial data can even be 
simulated without the damping at all. This is in sharp 
contrast to the typical values of the damping parameter, 
k ~ 1/A, required in simulations that employ the driver 
gauge conditions (|2"3"|) . (f2"Tj) . or 
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FIG. 17: A convergence in a 5D algebraic DWG evolution 
of $o = 0.4(1 — i), with h = 1/32. Top panels indicate 
that in the regions where a non-negligible amount of matter 
is present the convergence rate Q48p of the metric function A 
is essentially p ~ 2. Other functions, such as energy density, 
depicted in the bottom right panel, have similar convergence 
trends. The bottom left panel illustrates that the constraint 
violations consistently decrease with increasing resolution. 
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C. Convergence 



FIG. 18: The total residual (J49J) of the system of our FDA 
equations at the hypersurfaces z = 0, z = 1/2 and z = 1 for 
3 resolutions, with the coarsest one h — 1/32. The residual 
quickly decreases as a function of the numerical resolution. 



One of the crucial tests of numerical FDA schemes, 
such as one that we use here, involves the investigation 
of the convergence of the generated numerical solutions 



as a function of resolution. We perform convergence tests 
based on the assumption [43[ that for any of the unknown 
functions, Y(t,r,z), appearing in our system, the corre- 
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sponding discrete quantity, Yh(t, r, z) in the limit h — >• 
admits an asymptotic expansion of the form 

Y h (t, r, z) = Y{t, r, z) + h p e p {t, r, z) + ■ ■ • (47) 

where h is the spatial mesh size, e p (t,r,z) is an h- 
independent error function, and p is an integer that de- 
fines the order of convergence of the scheme. We consider 
sequences of three calculations performed with identical 
initial conditions, but with varying resolutions, h, h/2 
and h/A, and compute 

Myt^t)^' (48) 

for a grid function Y. 

We fix $ = 0.4 (1 — i) and use algebraic DWG evo- 
lution in 5D with the parameters = 3, q = 1/2, n = 
0.07, e = 0.125, and show in Fig. [T71 several plots illus- 
trating the convergence. The top left panel shows the 
radial dependence of the metric function A along z — 
at t ~ 2 M obtained in simulations using 3 different res- 
olutions, h,h/2 and h/A, with h = 1/32. The decreasing 
differences between solutions obtained using increasing 
resolutions indicate convergence. The top right panel 
in Fig. [T71 shows that the convergence rate (14"5)) is mostly 
quadratic, except in the asymptotic region where the am- 
plitude of the field is small and the computation is un- 
reliable. The bottom left panel depicts the Hamiltonian 
constraint at a late moment of the evolution: violations 
of the constraint are small and decrease further when 
the grid is refined. Finally, the bottom right panel illus- 
trates the time variation of the logarithm of the matter 
energy- density at (f, z) ~ (0.5,0). Again, a convergence 
is evident and is compatible with p^, 1.5 in the regions 
where a non- negligible amount of matter is present. 

An additional measure of accuracy of the scheme can 
be obtained by monitoring the total normalized residual 
of our equations that can be defined as 

{ n v n z j: y i ) l4yj 

where TZ-Yij is the residual of the FDA equation govern- 
ing the variable Y at a grid-point We show the 
logarithm of the residual in Fig. [18] for three resolutions 
h,h/2 and h/A, with h = 1/32 in the 6D simulation us- 
ing algebraic DWG and the initial scalar pulse amplitude 
$ = 1.0exp(-(r 2 + z 2 )/0.125). Evidently, the residual 
quickly decreases as a function of the numerical resolu- 
tion, again indicating convergence of the scheme. 

D. Conclusion 

We have described a generalized harmonic formula- 
tion of the Einstein equations in axial symmetry in 
D-dimensions and constructed the first numerical code 



based on it. We chose the coordinates in which the back- 
ground symmetries are explicit. This, however, resulted 
in a coordinate singularity on the axis, r = 0. While 
at the continuum level the equations of motion ensure 
regularity of a solution on the axis, extra care must be 
exercised so that this remains true in discrete numerical 
calculations. We have devised a regularization procedure 
that achieves that, while preserving the hyperbolicity of 
the evolution system. In our implementation we integrate 
the full D-dimensional equations and find that this ap- 
proach is smoother at the axis and is generically more 
stable compared to the approach that solves the 2+1 
equations, obtained by a dimensional reduction on the 
symmetry. 

We expect that our new code will enable systematic 
investigation of many problems of interest in axisymme- 
try. As a first application we tested the performance of 
the code in the context of fully nonlinear gravitational 
collapse of a complex, self-interacting scalar field propa- 
gating in a £>-dimensional Kaluza-Klein spacetime. We 
assumed spherical symmetry in the (D — l)-dimensional 
noncompact portion of the spacetime, which effectively 
reduced the problem to 2+1. The scenarios that we have 
considered ranged from the dispersion of dilute pulses to 
the collapse of strongly gravitating pulses that lead to 
black hole formation. One of the aspects of our code was 
the use of radial compactification which, in conjunction 
with sufficient Kreiss-Oliger-type dissipation, provided 
a viable alternative to the truncation of the spatial do- 
main and the use of approximate outer boundary condi- 
tions. Another ingredient of our algorithm that in some 
regimes was vital for long-term stability was the addition 
of constraint-damping terms to the evolution equations. 

We described several strategies to fixing the coordinate 
freedom that are compatible with the GH approach and 
experimented with those. Our studies of evolutions us- 
ing damped wave gauge that was enforced algebraically 
indicate robustness of this choice, its weak dependence 
on parameter settings, and stability. On the other hand, 
our experiments with various drivers — described in de- 
tail in Sec. IIII CI — reveal that in the case of symmetries 
these drivers are considerably less effective relative to the 
3 + 1 simulations that use Cartesian coordinates, a con- 
clusion similar to that drawn in (26j . However, it would 
be interesting to examine performance of the drivers in 
additional axisymmetric situations, other than collapse. 

One can consider several ways to improve the code. 
Specifically, it seems natural to use hyperboloidal slic- 
ing, similar to one suggested in [44j, instead of the spa- 
tial compactification that we presently employ. We ex- 
pect this will allow calculating asymptotic quantities and 
emitted gravity wave with a greater accuracy, and will 
help to improve late-time stability of the evolution, that 
is adversely affected by the unphysical radiation caused 
by the loss of numerical resolution near the outer bound- 
ary. One possible extension of our code which will 
broaden significantly the spectrum of possible axisym- 
metric configurations accessible with it, would be an ad- 
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dition of rotation. In addition, coupling the gravity dy- 
namics to evolution of more general matter, such as fluid, 
will enable studying certain situations relevant in astro- 
physics. 

Although the main purpose of this work was to de- 
scribe the code and test its performance in a non-trivial 
highly dynamical setting, our preliminary study of col- 
lapse in a Kaluza-Klein background indicates rich and 
distinctive phenomenology, deserving further investiga- 
tion. Among interesting questions, which will be ad- 
dressed elsewhere, is determining what classes of ini- 
tial data lead to formation of black holes of specific 
topology and constructing a phase diagram of the so- 
lutions (see [45| for some predictions concerning the dia- 
gram) . In addition, a detailed investigation of the situa- 
tion near threshold for black hole formation 12 and clas- 
sification of possible outcomes is necessary. In particu- 
lar, in this regime the fields typically oscillate, creating 
several outgoing shells. We found that the matter dis- 
tributes anisotropically within the shells, forming sepa- 
rate bounded systems, and that curvature singularities 
can develop inside those. We expect that the effect will 
be more pronounced in higher dimensions where many 
more shells can form near threshold [3lj . and since the 
gravitational field of an isolated system is increasingly 
localized in higher dimensions, the shells have a greater 
chance to separately become bounded systems capable of 
collapsing and forming black holes. It is worth exploring 
to which extent the expectations are true. However, such 
a study will require an improved horizon-finder that will 
be able to locate several moving horizons, and we are 
working in this direction. 

Finally, it would be very interesting to compute 
the detailed gravitational-wave signal emitted in var- 
ious regimes. While from the perspective of a four- 
dimensional observer all our solutions appear spherical 
that are not expected to emit any gravitational radia- 
tion, the waves are certainly produced and carry energy 
away. The "missing" energy, as it would be seen in 4D, 
will then provide a circumferential evidence in favour of 
the extra-dimensions. In addition, the spectrum of the 
gravitational waves must contain frequencies associated 
with the length scale of the extra dimensions. Therefore, 
by measuring the signal directly one can, in principle, 
probe the dimension and the topology of the spacetime. 
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Appendix A: Dimensionally reduced equations and 
boundary conditions 

Here we describe the approach that uses the Kaluza- 
Klein reduction of the D dimensional equations on the 
0(D — 2) symmetry and integrates the resulting three- 
dimensional equations. 

The most general D-dimensional line element that is 
invariant under action of the group of rotational symme- 
tries 0(D — 2) can be written as 



ds 2 



3 2aS dsl+e 2 f s dQl 



,2aS 0) 
9ab 



dx a dx b + e 2l3§ dnl. 



(Al) 



where the metric components are functions of three- 
dimensional coordinates x a alone, and a, j3 are constants, 
chosen for convenience below. 

We define the conformally rescaled metric g ab = 

e 2a ^g^ an d write the dimensional reduction of the 
Einstein-Hilbert Lagrangian as 

£eh = V~9dRd = 

R 3 + n{n-l)e- 2(l§ - 



-2n(3US -n{n + l) P 2 (d Sf 



where derivatives are computed using g 

as = 



(§sy = g^ ab d a Sd b S. 

Substituting the non-tilded metric and using the rela- 
tions between conformally related quantities [46| , 



{-h)- 1/2 da(g^ ab {-h) 1/2 d b s), 



(A2) 

(3)ab. 

and 



Ra = e 



as = 



dS) = 



-2aS 



(») 



-2aS 



-2aS 



R 3 -4anS-2a 2 (dSf 

a{dsf + as\ , 
(ds) 2 , 



(A3) 



where the derivatives in the right-hand-side are now com- 
puted with the non-tilded metric, we arrive at 



^-QdRd = V~93 e 



(a+ra/3)S 



R 3 - {Aa + 2n/3)nS 



(2 a 2 + 2 n a /3 + n(n + l)/3 2 ) (OS) 2 + 



n(n ~ 1) e 



-2^S+2aS 



(A4) 



By choosing a = —n /3, we convert the lower-dimensional 
action into Einstein-Hilbert form, C = \f— 53-R3. Since in 
this case DS is multiplied by a constant factor, it does not 
contribute to the equations of motion and can be omit- 
ted. In addition, fixing (3 — [n(n + l)] -1 / 2 ensures the 
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canonical normalization of the kinetic term of the scalar 
S. After these manipulations the Lagrangian becomes 



Ceh = V^9i R3 - (OS) 2 + n{n - 1) e 



-2c„5 



(A5) 



where we have defined c n = yin + This La- 

grangian describes three-dimensional gravity coupled to 
the self interacting 13 scalar field S. The components of 
the original _D-dimensional metric g^ (|A1|) are given in 
terms of the lower-dimensional fields as 



J D ) - P -(2/cn)S (3) 

y a b — e y a b ' 



9nn 



exp 



y/n(n + l) 



S 



9n- 



(A6) 



A reduction of the matter Lagrangian of our model ([2]) 
takes the form 



|a$| 2 + 2y(|$|)e- (2/c "^ 
yielding the total action of the system, 
5* = Seh + <S$ = 

"i? 3 - |d$| 2 - (dsf + 



(A7) 




+ n(n — 1) e 



-2c„S 



2V(|*|)e 



-(2/c„)S 



(A8) 



that describes two interacting scalar fields minimally cou- 
pled to gravity in three dimensions. Varying the action 
with respect to the 3-metric and the fields one obtains 
the equations of motion, 

Rab = T a b = T a b — g a0 T = 

= d {a $d b) $* + d a Sd b S- 
-g ab (-2 V e-W°«) s + n{n _ i )e ~^ n s 



as 



c n n(n- l)e 2c ™ 5 
dV 



+ — Ve 



(2/c„)S _ 



0, 



,-(2/ Cn )S = q 



(A9) 



Next, we bring the Einstein equations in this system into 
the GH form using the transformations analogous to ([5]) , 
([7]) but applied to the 3-metric, and add a constraint 
damping term similar to ©• The resulting hyperbolic 
system is evolved in time. 

We are interested in asymptotically Minkowski times 

a circle, M D_2a x S 1 solutions of (|A9j) . satisfying g^ -> 

i]ab, g^a -> r 2 , and $ -> 0. Using (jA6|) we find 

1 ^ Note that the potential term of the scalar S is proportional to 
the curvature of the n-sphere, and hence the scalar S is massless 



that asymptotic boundary conditions obeyed by the re- 
d uced field s in this case are g a 3 J — > ri ao r 2n and S — > 
y/n{n + 1) log(r). Since it is difficult to handle blowing- 
up conditions of this sort in numerical implementations, 
we redefine our variables by factoring out this singular 
behavior: g ao — > g a b"f 2n and S = S + y/n(n + 1) log(r). 
However, as a result of this transformation, the radial 
component of the GH source functions defined in (j4]) , ac- 
quires a term singular at the axis Hi = n/r + . . ., where 
ellipses designate regular at r — terms. This behavior is 
similar to what happens in the unreduced system and in 
analogy to that case we regularize the source functions 
by subtracting off this singular flat-background contri- 
bution, see Sec. IIII Al The gauge conditions are then 
applied to the regularized source functions. 

Regularity conditions at the center r = are again 
analogous to those in the unreduced case. Specifically 
3oo^ ' 9xi j 9o2 -i 522^ an d S are even functions in r as r — > 0, 
while and g\^ are odd. Moreover, requiring the ab- 
sence of conical singularity at r = places an additional 
condition that — exp(— c„ S) — 0(r 2 ) as r — »• 0. Since 
it is desirable to have a number of boundary conditions 
matching the number of dynamical variables, one might 
consider a regularization similar to (TT5)) that achieves this 
by defining A = (gn — e c ™ s )/r that behaves as A ~ 0(r) 
near r = 0, one than eliminates S from the scheme and 
uses A as a fundamental variable instead. However, a 
closer examination of the equation governing A reveals 
that the regularization does not work in this case be- 
cause the equation is not automatically regular at r = 
as it was in the unreduced approach. Rather the equation 
contains term proportional to 1 /r, which is regular only if 
the constraints are explicitly satisfied. Namely, not only 
the number of boundary conditions exceeds that of the 
dynamical fields — as it was before the regularization — 
but now the extra condition is also algebraically more 
complicated. Hence, instead of introducing A we opted to 
implement a more straightforward regularization method 
that maintains 5 as a fundamental dynamical variable 
and uses analytical Taylor-series expansion to compute 
its value near the axis; see 26] for more details. 

A comparison of the dimensionally reduced approach 
and the full D-dimensional method described in the 
main text, shows that both methods perform compara- 
bly in the weak-field regime when M ^ 0.5. However, for 
stronger initial data the dimensionally-reduced approach 
is considerably less stable and prone to developing co- 
ordinate singularities, regardless of the gauge conditions 
that we use. We do not fully understand the reason for 
this, a further investigation is required, and we hope to 
report on this puzzling issue elsewhere. 

in axisymmetry in 4D. 
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